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Mesoscopic quantum paradoxes, like the famous Schrodinger cat, are increasingly accessi- 
ble to experiment. While testing quantum mechanics under such conditions is an important 



chaUenge, the computation of ,uan,„ m predictions in these states is notorious* difncni. ow- 
ing to exponential complexity. The most straightforward approach of using probabilistic or 

3 ' 



'Monte-Carlo" sampling was thought by Feynman and others to be impossible, due to the fa- 



mous Bell inequality. It is an open question whether probabilistic simulations of these states 

in 

^ ■ can be carried out, and what are their error properties. Here we resolve this question by 
carrying out direct probabilistic simulations of several quantum paradoxes using a digital 
computer. We have treated multipartite Bell violations for up to 60 qubits, similar to those 
generated in photonic and ion-trap experiments, both with and without decoherence. Our 
results demonstrate that quantum paradoxes are directly accessible to probabilistic sampling 
methods, and we analyze the scaling properties of sampling errors. We anticipate that our 
results may provide a starting point for the development of more sophisticated computational 
tools, both for testing quantum theory at the boundary of the quantum and classical worlds, 
and for quantum engineering of new technologies that exploit this science. 
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The possibility of macroscopic entanglement^ of quantum states has led to many conceptual 
problems^. In order to understand this limit, Feynman asked: "Can quantum systems be proba- 
bilistically simulated by a classical computer?'^. His answer was "If. . . there's no hocus-pocus, 
the answer is certainly, No!". The argument rested on the key assumption that violations of phys- 
ically accessible Bell inequalities 4 could not be treated using probability distributions. Here we 
give counter-examples to this claim, through direct simulation of Bell's inequality violations. 

What is the technique that allows these probabilistic simulations of Bell violations? We sim- 
ulate quantum ensemble averages, not each individual observation. To achieve this, we employ 
the positive phase-space distributions of quantum optics^ for our calculations. These have statis- 
tical moments which correspond to those measured in Bell violations. We generate phase-space 
distributions for a number of quantum states used to investigate quantum entanglement and Bell 
inequalities. Our methods are applied to states of the type generated with both the photonic and 
ion-trap experiments. We demonstrate a clear violation of a Bell inequality in every case, showing 
the failure of any local hidden-variable theory. 

A Bell inequality is a constraint on observable correlations of a physical system that obeys 
a local hidden variable theory (LHV) 4 9 . This is a classical theory where measurements by two 
spatially separated observers (say, Alice and Bob) are obtained from random samples of a param- 
eter A. The measured values are functions of local detector settings and the hidden parameter A. 
Thus, the value observed by Alice with detector setting a is A(a, A), and similarly for Bob. All 
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correlations are obtained from a probabilistic calculation of form: 

C(X,Y)= [ X(X,a)Y(X,b)P(X)d\ (1) 

J A 

Here X(X, a) and Y(X, b) are experimental values, usually encoded as either 1 or —1 in a 
binary experiment. These assumptions lead to inequalities that any LHV correlations must satisfy. 
Quantum mechanics is known to violate these inequalities. But can one simulate these violations 
using probabilistic methods equivalent to quantum mechanics? Here we report probabilistic sim- 
ulations of Bell inequality violations for polarized photons (Fig. 1), as in the seminal experiments 
of Clauser 9 , Aspect et al.^and Zeilinger et alP-^, where: 

W = ^[|0>|1>|0>|1> + |1>|0>|1>|0>]. (2) 

To achieve this result, we use the positive P-representatiorP^, which is an exact expansion of 
an arbitrary density matrix of bosons with a positive probability distribution P(ct, (3). In terms of 
this general distribution, the correlation of quantum counts hi = ajSj at different locations is: 

where n £ = a^i. There is a remarkable similarity between the hidden variable theory © of Bell, 
and the positive-P formula © for quantum correlations. However, the positive-P representation 
can be used even for states that violate a Bell inequality^. The crucial reason for this difference is 
due to the different quantities sampled in the correlations. The fundamental observables in Bell's 
case, of form X(X, a), are equal to actual observed real numbers — i.e., (0, 1, . . .) for photon 
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counts. The corresponding observables in the positive-P case, of form n (a, (3), are complex num- 
bers whose mean values and correlations correspond to observable means and correlations. The 
results of the simulations are shown in Fig. 1, with the details in the Supplementary Methods on- 
line. These indicate a clear violation of a Bell inequality in both the two-particle case (N = 1) and 
the four-particle case (N = 2) which has also been observed experimentally^. 

To understand the ultimate scaling properties of probabilistic sampling methods, we have 
also simulated higher order correlations that violate multipartite Bell inequalities. These are found 
in quantum states that display Bell violations with M observers, not just two. The most well-known 
examples are the multimode entangled Greenberg-Horne-Zeilinger-Mermin (GHZM) states^L^, 
corresponding to "Schrodinger Cat" states. We therefore considered GHZM states which describe 
M spin-| particles or qubits: 

|$> = ^(|t...t> + e < *U--4»- (4) 

Here |t) and \\) denote spin-up or spin down particles in the ^-direction. As well as being 
of deep significance in quantum physics, such mesoscopic states have been generated in recent 
ion-trap experiments^^. Quantum paradoxes are obtained on measuring an operator A which is 
defined as a linear combination of 2 M distinct M -th order correlation functions: 

M 

For odd M we follow Merminpland take <j) = — 7r/2, Sj = 1 and 9j = 0. For even M we 
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follow ArdehalP^ and take <f) = ir, Sj = — 1 and 0j = for j < M, sm = L @m — — vr/4. Genuine 
multipartite Bell violations are known to exist for these states, which imply that the observed 
correlations cannot be explained by a nonlocality shared among M — 1 or fewer qubits. To test for 
genuine multipartite Bell nonlocality, we adopt the hybrid local-nonlocal LHV model introduced 
by Svetlichny^ and Collins et alM In order to simulate the Bell violation for these multipartite 
states, we use the SU(2)-Q distribution, which is well-suited to this type of Hilbert spaceP^. 

The difference between the positive-P and SU(2)-Q representation, compared to an LHV, can 
be illustrated by considering the case of M = 2, where simplifying the Ardehali inequality will 
give rise to the corresponding quantum mechanical expectation value: 

F^-i&l&D + i&l&l). (6) 

The correlation between real parts of the values of two factors for one of the terms is plotted 
in Fig. 2a and 2c, and the correlation between the two terms of © is plotted in Fig. 2b and 2d. 
Note that neither of the representations limits the values of Re ex* an d Re a 2 to the range [—1,1], 
as would happen in an LHV. This essential feature means that Bell's theorem does not limit our 
results, because the sampled values are not the same as their physical eigenvalues. This property 
of having a different apparent range from the operator eigenvalues is also shared by the theory of 
"weak" measurements^. Other examples of quantum correlations have been treated using these 
methods, including quantum solitons^, interacting quantum fields 25 equivalent to ~ 10 6 qubits and 
the Dicke superfluorescence modePS although no previous Bell violations were reported. These 
simulated correlations are in general agreement with experimental observations^. 
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The state © can be readily sampled using probabilistic random number generators together 
with the Q-f unction (details can be found in Supplementary Methods online). These samples can 
be used to calculate required expectation values (Fig. 3a). In the graph, the red dashed line is the 
minimum correlation required to demonstrate a Bell violation. Since the calculation is essentially 
a parallel one, we employed graphical processor unit (GPU) technology to calculate many samples 
in parallel, with a number of qubits ranging from M = 2 to M = 60. This corresponds to 
measurement of a quintillion (10 18 ) distinct sixtieth order correlation functions. Bell violations 
were verified in all cases, while genuine multipartite violations of LHV requiring all M observers 
to participate, with F sample / F QM > l/y/2, were verified for M < 50. 

To understand the source of sampling errors, we investigated the scaling of errors with 
system-size for single measurements of a low-order spin correlation (Fig. 3b), in addition to the 
exponentially large numbers of measurements in the expectation values F. As an example of 
low-order correlation we have chosen the total number of "spin-ups" iV = (J2jLi (H + 1) /2)- 
Low-order correlations were easily calculated with decreasing sampling errors as M increases. 
This explains why low order correlations could be efficiently sampled in previous work, for much 
larger Hilbert spaces than the ones studied here. 

In contrast to this, high-order correlations showed exponentially increasing sampling error. 
The relative error in F scales as 2 M / 3 , meaning that the time taken at constant error scales as 
22M/3 Surprisingly, probabilistic sampling scales more favorably than experiment, which would 
take time in proportion to 2 M , owing to the exponentially large number of measurements needed. 
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The increased efficiency of computational sampling is caused by the exponentially large number of 
measurements that are calculated in parallel in the phase-space simulations. This more than offsets 
the large number of samples required. 

Experimental measurements of such high-order correlations are limited due to noise and 
decoherence, which is an important issue in the observation of mesoscopic quantum effects 28 . The 
experimentally observed magnetic field noise found in ion-trap experiments can be easily added 
using a simple model. Following Monz et alP', we assume a delta-correlated noise such that 
(AE(t)AE(t')) = AEl5(t - t'), with an interaction Hamiltonian of the form: 

By analyzing operator correspondences in the Heisenberg picture, we find that this can also be 
readily simulated. This was achieved by multiplying an independent noise term exp (ieMQ) by 
the value corresponding to the operator A in each of the samples after every time step At. Here 
e = AE \/At/h defines the speed of the decoherence, and Q is a Gaussian random number such 
that ((j(j') = Sjji . This is shown in Fig. 4, which demonstrates the experimentally observed 
quadratic super-decoherence as M increases. 

We see from the present results that while higher order correlations in GHZM states can be 
simulated, they require many samples to reduce errors to acceptable levels. This result is not com- 
pletely unexpected. Higher-order correlations are often associated with isolated quantum states, 
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which effectively become "needles in a hay-stack" in large Hilbert spaces. Universal digital quan- 
tum computers are expected to become useful to overcome this problempl. However, these are 
currently limited in size to around 6 qubits or less^l. Scientifically, to test quantum theory it is also 
important to obtain quantum predictions using methods that do not inherently rely on quantum 
mechanics being correct. 

In summary, it is not impossible to use probabilistic digital algorithms to simulate meso- 
scopic quantum paradoxes comparable to those found in current experiments. Such quantum tech- 
nologies are already finding applications in quantum logic based encryption, atomic clocks and 
interferometers. Digital simulations could have a direct application to the design and implemen- 
tation of these technologies, as well as providing quantum predictions for fundamental tests of 
quantum physics. 
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Supplementary Information 

* Supplementary Methods 

In this supplementary information, we give the technical details of the probability distribu- 
tions and sampling techniques used to obtain the data in the main paper. 
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Bell inequalities. The sampling of the Bell inequalities is performed using the positive-P repre- 
sentation of the correlated multi-photon target state. Here we consider correlated states with N 
photon pairs, which are described by the following ideal quantum state: 



(a[ + a\ + + a\_a\\ |0) 

\N) = ± , (8) 

iV!(jV + l) 1/2 

where N is the number of photon pairs. These states are generated in certain types of optical 
parametric down-conversion experiments, and have been studied by several different groups. We 
define the intensity correlations to be G /J (7, N) = (iV|y4 7 (l, ai)A J {^/, a 2 )\N), where 7 is a linear 
combination parameter related to a polarizer angle, and auxiliary functions are introduced so that: 



i J ( 7 , o fc ) = (Vt4- + V 1 ~ t4+) J (VTOfe- + \fi - 7flfc+) 3 , (9) 



i J (oo, a k ) = : (a\_a k _ + a\ + a k ^j :. (10) 

Tthe standard two-particle CHSH form of Bell's inequality given by Clauser, Home, Shi- 
mony and Holt, is especially important, as it gives classical limits to the expected correlation for 
the above experiment: 

C[A(a),B(b)} + C[A(a),B(b')} + C[A(a'),B(b)] - C[A(a% B(b')} < 2. (11) 

More generally, for iV photon pairs, the Bell-type inequality is then known to be of the form: 
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A(9) = 3g(9)-g(39)-2<0, 



(12) 



where 

g J N (B) = G JJ (cos 2 9, N)/G JJ {oo, N). (13) 

While this expression generalizes the usual CHSH and Bell expressions to a multi-particle form, 
the quantum mechanical prediction for g J N has especially simple form of g^{9) = cos 2N 9 for the 
cases J = N, N = 1,2 which we have simulated. This gives violations of 

Ai pair (0) = 3cos 2 0- cos 2 39 - 2, (14) 

for the two particle case, which corresponds to the usual two-particle experiment originally pro- 
posed by Bell. For the four-particle case one finds that: 

A 2pairs (#) = 3 cos 4 0- cos 4 30 -2. (15) 

In all cases, the corresponding photonic system is a four-mode state \N), which has the 
corresponding positive-P distribution 
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P(a,(3) 



(P 1+ + a* 1+ ) (&+ + a* 2+ ) + (A- + <j (fe- + tg_) | 2jV 
(2vr) 8 (iV + l) (iV!) 2 2 4iV 



exp 



(16) 



This distribution exists and has a positive, probabilistic behavior for all values of N. We sampled 
this distribution by transforming to a form where we could use the well-known von Neumann re- 
jection algorithm, combined with an underlying lambda-distribution. The exact results are plotted 
in Fig. 1, together with the probabilistic sampling results. 



GHZM states and SU(2)-Q representation To handle systems with many distinct observers — 
known as multipartite systems — we found it most efficient to use a representation that was op- 
timized for treating quantum binary states or qubits. These are the SU(2) coherent states and the 
corresponding positive Q-function. For an M-mode system the un-normalized SU(2) coherent 
states are defined as: 

M 

(17) 



l»> = e 5t -|0) = n |0>i + ^|l>j 



The SU(2) Q-function is the expectation value of the density operator over the SU(2) coher- 



ent states, and is defined explicitly in a normalized form as: 



' M 

n 



7T 



2\3 



(i + N 2 ) 



( z IIpII z ) 



(18) 



For the GHZ states ©, the SU(2) Q-function is easily computed using the orthogonality of the 
spin basis. This leads to a positive distribution for a GHZ state with M sites: 



1 M 



2\3 



2 j=1 TT (1 + \ Zj \ 2 ) 
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A I 



Hzj + e 



(19) 



In this case the sampling was again carried out by using a von Neumann rejection algorithm com- 
bined with an inverse cumulative distribution which was sampled exactly. 

The expectation value of the operator A © which gives the multi-partite Bell violations 
in terms of the SU(2)-Q function is also able to be computed, using the properties of the spin 
operators. It is given by: 

f ™ 3 ((1 + Sj )z* + (1 - Sj)zj) e~ is ^ 
(<S>\A\<t>) = /dzQ(z)n— 3) '- \ J 3) • (20) 

J fJl l + \Zj\ 

The expectation of the number of "spin-ups" can be shown to be 



.7=1 J .7=1 V 
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1 + 


Vi 


2 



(21) 



As described in the main text, for odd M we follow Mermin 15 and take 



-7T 



1% Sj 



and 6j = 0. The value of interest in this case is F = Im($|A|$), which is limited to 2( M-1 " 2 
by any LHV, but has a quantum mechanical (QM) expectation of 2 M_1 . For even M we follow 
Ardehali^and take <p = tt, Sj = —1 and 9j = for j < M, sm = 1, &m = — 7r/4. The value 



of interest is F 



Re($|i|$) + Im($|i|$) with a QM expectation of 2 M ~ 1 / 2 . Here the LHV 



limit is 2 A/ / 2 . Therefore, for any M we have a relative violation 2^ M 1 ^ 2 . 

This result is used to calculate the variance in the low-order correlations. All of these 
results were calculated using a numerical code that generates the samples, calculates the cor- 
responding correlation function, and accumulates the result for averaging. In order to obtain 
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low enough sampling error to give the strong correlations needed to violate 60th order Bell in- 
equalities, and show genuine multipartite Bell violations, we needed to use adequate numbers of 
samples. The results given here required 10 12 samples, which were obtained using GPU multi- 
processor technology. The results reported here required up to about 24 hours on four paral- 
lel GPU processors, each with 448 cores. The code for the simulations is available publicly at 
http : / / git hub . com/Mant icore /bells im- letter. 
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Figure 1 Bell violation. Simulated Bell violation A as a function of the relative polarizer 
angle 9 for one (a) and two (b) photon pairs using the positive-P distribution with 2 21 (a) 
and 2 25 (b) samples. The filled area corresponds to the estimated error range around the 
mean of A(6>) for the sampled state, while the dashed line is the exact quantum mechani- 
cal prediction of this value. 
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Figure 2 Spin moment correlations. Correlations for the different parts of the quan- 
tity (HI), in case of positive-P (a, b) and SU(2)-Q (c, d) representations, 10 8 samples. 
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Figure 3 Violations for multi-particle GHZ states, a, Simulated Mermin violation us- 
ing SU(2)-Q representation. The values of expectations and errors are normalized by 
the quantum mechanical prediction for the corresponding M. The horizontal grey dashed 
line gives the quantum prediction. The error bars show the sampled result and estimated 
sampling errors at each value of M. The red dash-dotted line is the LHV prediction, which 
gives a Bell violation when above this line. Genuine multipartite Bell violations occur for 
-^sampic/^QM > 1/V2. b, Relative errors for F (blue line) and first order correlation, or 
total number of "spin-ups" (green dashed line) using SU(2)-Q representation. The dotted 
reference line shows the point at which the sampling errors would give scaling properties 
as slow as an experimental measurement. 
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Figure 4 The effect of decoherence on measured quantities. Decay of the sam- 
pled quantity F after the application of the simple model of super-decoherence (0), for 2 
(solid blue line), 3 (red dashed line), 4 (green dash-dotted line) and 6 (yellow dotted line) 
particles, with decoherence rate e = 0.1. The horizontal axis is the dimensionless time, 

r = t/At. 
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